Similarly to Notebook 2, we want to learn about the optical properties of the resonator that is causing the fringes we see in MRS spectra. Recall that in Notebook 2 the resonator was a Fabry-Pérot etalon which was (safely) approximated as one element with two plane-parallel mirrors, made out of the same material. Now the resonator under scrutiny is the MIRI short-wavelength detector. Contrary to a Fabry-Pérot etalon, the MIRI detectors are made out of a number of consecutive optical layers, each with their own geometric thickness and optical properties (more about that in Notebook 8). But, as such, it is no longer safe (i.e. justifiable) to model the fringes using a single reflectance Fabry-Pérot transmittance function. In this notebook, we do so anyway. The reason is, unjustifiable as it may be, such an analysis can yield zeroth order constraints on the optical properties of the detector (a.k.a. we need to start from somewhere).
import funcs
import mrsobs
import numpy as np
from astropy.io import fits
from scipy.optimize import curve_fit
from scipy.interpolate import InterpolatedUnivariateSpline
from matplotlib import pyplot as plt
%matplotlib notebook
import warnings
warnings.simplefilter('ignore')
We load the images for one band of the MRS for different kinds of sources, including:
Additionally the pixel-to-wavelength calibration map and the pixel-to-along-slice position map are imported.
# Define paths to data
workDir = '/Users/ioannisa/Desktop/python/miri_devel/'
cdpDir = workDir+'cdp_data/'
d2cMapDir = workDir+'distortionMaps/'
lvl2path = workDir+'FM_data/LVL2/'
# Get data
band = '1A'
ext_source_sci,ext_source_bkg = mrsobs.FM_MTS_BB_extended_source(lvl2path,band,bb_temp='800K')
ext_etal_source_sci,ext_etal_source_bkg = mrsobs.FM_MTS_800K_BB_extended_source_through_etalon(lvl2path,band,etalon='ET1A')
point_source_sci_p1,point_source_bkg_p1 = mrsobs.FM_MTS_800K_BB_point_source_raster(lvl2path,position='middle',pointing='P1')
# Get wavelength calibration pixel map
d2cMaps = funcs.load_obj('d2cMaps_band{}'.format(band),path=d2cMapDir)
lambdaMap = d2cMaps['lambdaMap']
lambdaMap[lambdaMap==0] = np.nan
alphaMap = d2cMaps['alphaMap']
nslices = d2cMaps['nslices']
We subtract background exposures where available
# perform transform
ext_source_bkgsubtr = ext_source_sci-ext_source_bkg
ext_etal_source_bkgsubtr = ext_etal_source_sci-ext_etal_source_bkg
point_source_p1_bkgsubtr = point_source_sci_p1-point_source_bkg_p1
We perform an even-odd row signal correction to the data (caused by the read-out pattern of MIRI detector pixel rows)
ext_source_oddevencorr = funcs.OddEvenRowSignalCorrection(ext_source_bkgsubtr)
ext_etal_source_oddevencorr = funcs.OddEvenRowSignalCorrection(ext_etal_source_bkgsubtr)
point_source_p1_oddevencorr = funcs.OddEvenRowSignalCorrection(point_source_p1_bkgsubtr)
We perform the same analysis as in Notebook 2, but rather than carrying out the analysis for etalon data we carry it out for extended source fringes. Although the latter have much smaller contrast compared to the former (i.e. we don't see sharp etalon lines but a blend of etalon lines), the working concepts are exactly the same. So let's try it, let's start by fitting a spectral profile to the fringe peaks to compensate for the sampling effects.
# Pixel trace in MRS slice
ypos,xpos = funcs.detpixel_trace(band,d2cMaps,sliceID=nslices/2,alpha_pos=0.)
# Extended spectrum fringe peaks
ext_source_peaks = funcs.find_peaks(ext_source_oddevencorr[ypos,xpos], thres=0.3, min_dist=6)
# Fit etalon lines
fitparams,errors,fitting_flag,range_ini,range_fin = funcs.fit_etalon_lines(lambdaMap[ypos,xpos],ext_source_oddevencorr[ypos,xpos],ext_source_peaks,fit_func='gauss1d',sigma0=0.01)
lineheights = funcs.get_amplitude(fitparams=fitparams,fitting_flag=fitting_flag)
linecenters = funcs.get_linecenter(fitparams=fitparams,fitting_flag=fitting_flag)
# Take univariate spline through the fitted peaks of the etalon lines
interpolator = InterpolatedUnivariateSpline(linecenters,lineheights,k=3,ext=3)
ext_source_peakprofile = interpolator(lambdaMap[ypos,xpos])
# let's look at the result
fig,axs = plt.subplots(3,1,figsize=(12,13))
for plot in range(2):
axs[plot].plot(lambdaMap[ypos,xpos],ext_source_oddevencorr[ypos,xpos])
for i in range(len(fitparams)):
plotx = np.linspace(fitparams[i][1]-10*fitparams[i][2],fitparams[i][1]+10*fitparams[i][2],100)
ploty = funcs.gauss1d_woBaseline(plotx,*fitparams[i])
axs[plot].plot(plotx,ploty,'r')
axs[plot].plot(linecenters,lineheights,'ro')
axs[plot].plot(lambdaMap[ypos,xpos],ext_source_peakprofile)
axs[plot].set_ylabel('Signal [DN/sec]',fontsize=20)
axs[2].plot(lambdaMap[ypos,xpos],ext_source_oddevencorr[ypos,xpos]/ext_source_peakprofile)
axs[1].set_xlim(5.18,5.43)
axs[2].hlines(1,4.88,5.77,'orange')
axs[2].set_ylim(0.85)
axs[2].set_ylabel('Normalized signal',fontsize=20)
for plot in range(3):
axs[plot].set_xlabel('Wavelength [micron]',fontsize=20)
axs[plot].tick_params(axis='both',labelsize=20)
plt.tight_layout()
Let's perform the same exercise, but this time use the fringe pixel (/sampled) peak as the nominal fringe peak. This will of course result in sustaining aliasing effects in the final fringe transmission spectrum. However we base our actions on the conclusion that fitting a spectral profile to the fringe peaks does not yield better results.
# Pixel trace in MRS slice
ypos,xpos = funcs.detpixel_trace(band,d2cMaps,sliceID=nslices/2,alpha_pos=0.)
# Normalize the extended source spectrum to the fringe peak profile
ext_source_norm = funcs.norm_fringe(ext_source_oddevencorr[ypos,xpos], thres=0.3, min_dist=6)
# let's look at the result
fig,axs = plt.subplots(2,1,figsize=(12,8))
axs[0].plot(lambdaMap[ypos,xpos],ext_source_norm[0])
axs[0].plot(lambdaMap[ypos,xpos][ext_source_norm[1]],ext_source_norm[0][ext_source_norm[1]],'ro')
axs[0].plot(lambdaMap[ypos,xpos],ext_source_norm[2])
axs[1].plot(lambdaMap[ypos,xpos],ext_source_norm[0]/ext_source_norm[2])
axs[1].hlines(1,4.88,5.77,linestyle='dashed')
axs[0].set_ylim(21)
axs[1].set_ylim(0.85)
axs[0].set_ylabel('Signal [DN/sec]',fontsize=20)
axs[1].set_ylabel('Normalized signal',fontsize=20)
for plot in range(2):
axs[plot].set_xlabel('Wavelength [micron]',fontsize=20)
axs[plot].tick_params(axis='both',labelsize=20)
plt.tight_layout()
Let's carry on, as per Notebook 2, with Test 1.
Information about the optical properties of the detector can be attained by studying the fringes in wavenumber space. [Recall:] Assuming that the detector layers are plane-parallel, and that the incidence angle is small, subsequent peaks in transmission occur at a constant distance from one another. This distance relates to the optical thickness of the resonating layer.
fringepeaks_wavelength = lambdaMap[ypos,xpos][ext_source_norm[1]] # microns
fringepeaks_wavenumber = np.flipud(10000./fringepeaks_wavelength) # cm-1
mean_fringepeak_separation = np.mean(np.diff(fringepeaks_wavenumber)[1:-1]) # omit first and last data point
# fit straight line through distance data for comparison
popt,pcov = curve_fit(funcs.straight_line,fringepeaks_wavenumber[1:-2],np.diff(fringepeaks_wavenumber)[1:-1])
print r'Mean etalon line separation in wavenumber space is: Δσ = {} cm-1'.format(round(mean_fringepeak_separation,2) )
print r'The optical thickness D of the etalon is related to Δσ as: D = 1/(2Δσ) = {} um'.format(round(10000./(2*mean_fringepeak_separation),2) )
# Let's look at the results
fig = plt.figure(figsize=(12,5))
axs1 = fig.add_subplot(111)
axs2 = axs1.twiny()
axs1.plot(fringepeaks_wavenumber[1:-2],np.diff(fringepeaks_wavenumber)[1:-1])
axs1.plot(fringepeaks_wavenumber[1:-2],funcs.straight_line(fringepeaks_wavenumber[1:-2],*popt))
axs1.hlines(mean_fringepeak_separation,fringepeaks_wavenumber[0],fringepeaks_wavenumber[-1],linestyle='dashed')
axs1.set_xlabel(r'Wavenumber $\sigma$ [cm-1]',fontsize=20)
axs1.set_ylabel(r'$\Delta \sigma$ [cm-1]',fontsize=20)
axs1.tick_params(axis='both',labelsize=20)
axs2.set_xlim(axs1.get_xlim())
tickmarks = np.array([1750.,1800.,1850.,1900.,1950.,2000.,2050.])
axs2.set_xticks(tickmarks)
axs2.set_xticklabels(funcs.tick_function(tickmarks))
axs2.set_xlabel('Wavelength [micron]',fontsize=20)
axs2.tick_params(axis='both',labelsize=20)
plt.tight_layout()
Before, for the etalon data, we assumed that the etalon's reflecting mirrors were made out of the same material, and thus a single (wavelength-dependent) intensity reflectivity could be determined for the optical element. Although we know this to not be true for the detector (zinc-sulfide for the anti-reflection coating and aluminium for the pixels), we can still determine an "effective" reflectivity from the measurements. This can simply be done by looking at the fringe transmission trough level and relating that to the etalon finesse value, thus reflectivity.
# at lower wavelength end
T_e1 = 0.933
F,R_solution1,R_solution2 = funcs.reflectivity_from_continuum(T_e1)
print ' Finesse = {} \n R1 = {} \n R2 = {} (physically valid solution)\n'.format(F,round(R_solution1,2),round(R_solution2,2) )
R_low = R_solution2
# at higher wavelength end
T_e2 = 0.882
F,R_solution1,R_solution2 = funcs.reflectivity_from_continuum(T_e2)
print ' Finesse = {} \n R1 = {} \n R2 = {} (physically valid solution)'.format(F,round(R_solution1,2),round(R_solution2,2) )
R_high = R_solution2
plt.figure(figsize=(12,4))
plt.plot(lambdaMap[ypos,xpos],ext_source_norm[0]/ext_source_norm[2])
plt.hlines(1,4.88,5.77,linestyle='dashed')
plt.hlines([T_e1,T_e2],4.88,5.77,'r',linestyle='dashed')
plt.ylim(0.85)
plt.xlabel('Wavelength [micron]',fontsize=20)
plt.ylabel('Normalized signal',fontsize=20)
plt.tick_params(axis='both',labelsize=20)
plt.tight_layout()
We use the information derived from the extended source data to model the transmittance function of the detector. It is important to note which dataset the modeling is based on. We know that the fringe transmission differs with different optical stimuli and their spatial properties!
wavelengths = np.linspace(4.88,5.77,10000) # um
wavenumbers = np.flipud(10000./wavelengths) # cm-1
R = np.flipud(np.linspace(R_low,R_high,10000)) # [-]
D = 1./(2*mean_fringepeak_separation) # cm
phi = np.pi # rad
T_e = funcs.FPfunc(wavenumbers,R,D,phi) # [-]
# Let's look at the result
fig,axs = plt.subplots(2,1,figsize=(12,8))
for plot in range(2):
axs[plot].plot(lambdaMap[ypos,xpos],ext_source_norm[0]/ext_source_norm[2])
axs[plot].plot(wavelengths,np.flipud(T_e))
axs[plot].set_ylim(0.85)
axs[plot].set_xlabel('Wavelength [micron]',fontsize=20)
axs[plot].set_ylabel('Normalized signal',fontsize=20)
axs[plot].tick_params(axis='both',labelsize=20)
axs[1].set_xlim(5.503,5.692)
axs[1].set_title('<Zoomed view>',fontsize=20)
plt.tight_layout()
Here we attempt to fit the fringes with an optimized transmittance function. We do this by scanning through the fringe transmission profile, with a scan width equal to the period of the fringes.